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ABSTRACT 

Context. When inverting solar spectra, image degradation effects that are present in the data are usually approximated or not consid- 
ered. 

Aims. We develop a data reduction method that takes these issues into account and minimizes the resulting errors. 

Methods. By accounting for the diffraction PSF of the telescope during the inversions, we can produce a self-consistent solution that 

best fits the observed data, while simultaneously requiring fewer free parameters than conventional approaches. 

Results. Simulations using realistic MHD data indicate that the method is stable for all resolutions, including those with pixel scales 

well beyond those that can be resolved with a 0.5m telescope, such as the Hinode SOT. Application of the presented method to reduce 

full Stokes data from the Hinode spectro-polarimeter results in dramatically increased image contrast and an increase in the resolution 

of the data to the diffraction limit of the telescope in almost all Stokes and fit parameters. The resulting data allow for detecting and 

interpreting solar features that have so far only been observed with Im class ground-based telescopes. 

Conclusions. A new inversion method was developed that allows for accurate fitting of solar spectro-polarimetric imaging data over 
a large field of view, while simultaneously improving the noise statistics and spatial resolution of the results significantly. 

Key words. Techniques: imaging spectroscopy, polarimetric, methods: data analysis, numerical 



1. Introduction 

In solarphysical research, as in the vast majority of astrophysical 
research, information about the object under investigation is re- 
ceived in the form of electromagnetic radiation. Since this radia- 
tion was emitted by, passed through or reflected ofl' the object of 
interest, interpreting the information contained in the radiation 
typically involves computing the observed data from a model of 
the object using basic physical principles. 

Interpreting the data typically starts by making an educated 
guess about the observed object's atmosphere (the model), fol- 
lowed by an evaluation of the compatibility of the observable, 
calculated using this model, with the observed data. Additional 
improvements may then be made, depending on the flexibil- 
ity of the model. Optimizing the compatibility of the syn- 
thesized results with the data is usually conveniently formu- 
lated in terms of the minimization of the square of the dif- 
ferences. A number of codes exist that automate this opti- 
mizati on process, either by a standard lin earized opti mization 
(SIR (lRuizCo bo&delToroIniestal ll992h . SPINOR (iFrutigei 



20001: iFrutiger et al.l l2000h. N jCOL^ JSocas-Navarro et al 



19981) . HAZEL (! Asensio Ramos e t al."2008\, and many more); 
by statistical methods (HELIX (Laggetal. 2004, 2009), 
MERLIN, etc.); by means of Bayesian methods (Bayes-ME 



lAsensiq Ramos et aP l2007l): or by using artifici al neural net- 
works (ISocas-Navarrol 120031: ICarroU et al]l2008h . with greatly 
varying levels of detail in their treatment of the physics and often 
with the aim of treating very specific types of data. 

Usually the simplest atmospheric model that can reproduce 
the observations is used for analysis of the solar lower atmo- 
sphere. The model is characterized by a number of nodes, at 



which the physical parameters that are relevant for the forma- 
tion of the spectrum are defined. These quantities are then inter- 
polated to a finer grid to allow for the calculation of an accurate 
solution to the radiative transfer equation. Automatic inversion 
codes, which optimize this kind of model, have been in use for 
two decades for interpreting solar data. 

Spurred by the increase in computing power in recent times, 
much efifort has gone into increasing the level of detail and 
self-consistency in the treatment of the basic physics, non-LTE 
(NICOLE) and scattering polarization eflTects (HELIX, HAZEL, 
etc.). Without these, many observed spectral features can be 
shown to be poorly reproduced, even in cases where the fitted 
atmosphere is realistic. 

The observed data used in this process, however, quantify 
the properties of radiation that has passed through a telescope, 
an instrument, and a detector, all of which alter the data depend- 
ing on the size, optical quality and other instrumental properties. 
For ground-based observations, in addition, the radiation passes 
through the earths atmosphere, which redistributes it by refract- 
ing and scattering. 

A large body of work exists on removing the degrada- 
tion of the spatial resolution due to an increase in the ex- 
tent of the point spread function from instrumental and at- 
mospheric efifects. In particular the optimization of image in- 
formation from observe d data degraded by a tmospheric tur- 
bulence ((M O)MFBD (Ivan Noort et al.1 12005]), phase diver- 
sity methods (iPaxman et all 1 19921: lL5fdahl & Scharme3 11994 
and references therein), and speckle reconstruction methods 
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(Labevrie 1970"; "Knox & Thompson 1974; iLohmann etaDll983l: 
deBoeretal. 1992: von derLuehe 1991 and many more)) 
has recently found widespread application in spectral imaging 
data produced by an array of narrow band filter instruments; 
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e.g. S OUP (iTitle & Rosenberdll981h . CRISP^harmer 2006), 
GFPI (iPuschmann et al.1 l2006l) and IBIS (ICavallinii i2006b . to 
name just a few. These methods only work well for variable 
degradation on timescales [much] shorter than the time resolu- 
tion of the observation. 

In the spectral dimension, on the other hand, degradation by 
instrumental and atmospheric effects has a rather different effect 
on the data. Where the spatial degradation merely decreases the 
resolution in the spatial dimensions, the blending of the spec- 
tra results in observed spectra containing a mix of information 
from a collection of different atmospheric components. Much 
work has been done on the analysis and in terpretati on of such 
''multi-component" spectra (^Stenflol 1 19681: ISkumanich & LitesI 
ernascom & SolankilTT996r and many more), by mod- 
eling them as arising from two or more separate atmospheric 
''c omponents", which remain un resolved in the data. 

lOrozco Suarez et aP (l2007l) model the effect of telescope 
diffraction on the spectra as a local "stray-light" contribution, 
where the light from the surroundings is taking the place of the 
contribution of a second atmospheric component. The apparent 
multi-component spectra are not interpreted in terms of unre- 
solved fine structure, but rather in terms of an additive contribu- 
tion produced by telescope diffraction. 

In this paper we develop this idea further and develop a 
method that rigorously deals with the effect of image degrada- 
tion by instrumental and telescope effects, by integrating these 
effects in a 2-D adaptation of the SPINOR inversion code. 

2. Inversion of spatially degraded data 

Any process that infers the most probable atmospheric structure 
that gave rise to the observed spectral data may generally be re- 
ferred to as an inversion process. One could even argue that the 
inversion is only defined by the way the probability of a partic- 
ular atmospheric structure is defined, since the method used to 
locate the point of maximum likelihood is normally quite inde- 
pendent of the formulation of the probability function, but this 
would not do justice to the information that must be obtained 
from the atmospheric structure in order to optimize it. 

Considerable differences exist between the various inversion 
codes in the way they describe the atmosphere, the description 
of the radiative transfer problem, and in the way they try to max- 
imize the likelihood function. However, what they all have in 
common is that ultimately they all compute synthetic profiles 
from estimated atmospheres and compare these to the observed 
profile, in order to infer the most likely atmosphere that pro- 
duced the observed spectrum. 

A physically accurate description is included in most of these 
codes in ample detail, since it is obvious that this is needed to 
have any chance of interpreting an observed spectrum in terms 
of physical quantities. Although the most common assumption 
made in inversion codes is that of local thermodynamic equi- 
librium (LTE), there are codes that compute the line profiles 
in full NLTE (e.g. NICOLE, HAZEL), codes that include de- 
tailed scattering physics (HELIX, HAZEL), and codes that in- 
clude detailed molecular physics (SPINOR). In addition, there 
is a great variation in the amount of detail in which the atmo- 
sphere is parameterized, from slabs with constant properties to 
full 3-D structures with a full height stratification in all ph ys- 
ical quantities (see for instance lAsensio Ramos et aL 1 120T2L for 
an overview). 

Instrumental effects, on the other hand, with the exception of 
spectral smearing, are not generally considered in much detail, 
if at all, by many codes, because this is considered by most code 



authors to be taken care of by appropriate data reduction rou- 
tines. Although careful calibration of the data generally receives 
a lot of attention, there are always significant limits to the extent 
to which it is possible to remove the instrumental effects from 
the data. 

In particular, the removal from the data of any degradation 
of the data that reduces the information content (such as con- 
volution with an extended convolution kernel) will result in am- 
plification of any additional modification that was made to the 
data after that degradation took place, regardless of whether this 
modification is in the form of detector or photon noise, detector 
nonlinearities, or inhomogeneities in the response of any of the 
elements taking part in the acquisition of the data (such as filters 
or modulators). For this reason, it is advantageous to apply the 
instrumental effects to the synthetic data as a part of the relevant 
physics, rather than attempting to remove them from the data it- 
self. The most common instrumental degradation treated in this 
way so far is the degradation in the spectral dimension. Since 
the effect of this degradation on the spectra is frequently very 
severe, convolution of the synthesized spectra with the instru- 
mental transmission profile before comparing to the observed 
data is a standard ingredient in almost any inversion code. 

A rudimentary treatment of spatial degradation effects can 
be found in the form of a micro- and/or macro-turbulent broad- 
ening of the spectral line profiles in almost all inversion codes; 
however, a thorough treatment of degradations in the spatial do- 
main has thus far not been employed. This is probably not due 
to the relative importance of spatial degradation as compared to 
spectral degradation, but more likely due to the limitations im- 
posed by the available datasets. Recent advances in instrumen- 
tation (Hinode SP, IBIS, CRISP, etc), have produced datasets 
sampling both the spatial and the spectral dimensions with high 
enough resolution to open up the possibility of improved treat- 
ment 

lOrozc o Suare z et aO (l2007l) used the observed data from 
neighboring pixels to assemble a "local stray-light" spectrum, 
which contributes to the observed spectrum in a particular pixel 
of Hinode SP slit-scanned data. As acknowledged by the authors, 
this method is far from exact, but is simple to understand and 
easy to implement, and it fixes most of the fitting problems they 
encountered. 

However, since the estimate of the "stray-light" contribution 
can only be calculated from the observed data, this method in ef- 
fect approximates a deconvolution by subtracting a running aver- 
age. Apart from being formally incorrect, this procedure has the 
notable disadvantage that it removes an estimate of the signal not 
originating in a particular pixel (local stray-light) from the ob- 
servations, thereby substantially decreasing the signal-to-noise 
ratio of the data. In addition, the process of using the data to cal- 
culate a fit parameter (the stray- light fraction) makes the solution 
vulne rable to systematic errors (lAsensio Ramos & Manso Sainzl 

[mil). 

In the following, we develop an integrated approach to the 
problem of inverting spatially and spectrally degraded data. To 
see how to do this, we must first analyze the problem in some 
more detail. 

2.1. Image response 

The collection of solar spectral imaging data is no different from 
any other imaging problem, in that an image of the solar surface 
is formed by a telescope through an aperture in a focal plane, 
where it is sampled by an instrument that can detect both the 
spatial and spectral variability of the incident radiation. As a 
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consequence of the finite aperture of any imaging device, the ra- 
diation originating in a point on the object is redistributed across 
a spatially extended region on the detector. 

The implication of this for two neighboring points of the ob- 
ject is that their respective spectra have been averaged before 
they are detected in the focal plane. In the probable case that the 
spectra diff'er, the result is a mixed spectrum that may not even 
resemble the original spectra. 

To describe the eff'ect of the spatial smearing with a PSF 
(f)(x,y,A), we consider the data observed if only a single point 
source is present in (xq, yo)'- 

d(x, y, A) = i(xo, yo, A)(p(x - xo,y -yo,A). 

For a spectrum that is completely described by a set of parame- 
ters Of/ we may then write for the response of the observed data 
to a change of the parameter of/ in location (x, y)\ 

dd(x,y,A) dI(xo,yo,A) 

= (p(x- xo,y-yo,A) (1) 

oai oai 

since the instrumental PSF does not depend on the at. 
We recognize the derivative on the RHS as the response 
function of the modeled spectrum used in most inversion 
methods employin g a "greedy" algorithm (see for instance 
ICornien et al.ll200ll for a description) to minimize a merit func- 
tion (iLandi Degl'Innocenti & Landi Degl'Innocentill 19771) . The 
above expression gives us the true response of the data, given 
the undegraded response of the spectrum. We use this below to 
formulate the correct response matrix for a 2-D field of point 
sources, degraded by a known PSF. 

Although not strictly necessary, in the following we assume 
that the spatial part of the PSF does not depend on the wave- 
length so that we can write 

A) = (p{x,y)il/{A). 

There are many situations where this assumption is not valid, 
but it allows us to carry out the convolution over the wavelength 
before considering the spatial coupling, which makes the imple- 
mentation somewhat simpler. We intend to relax this restriction 
in future work. 

2.2. Levenberg-Marquardt minimization 

The Levenberg-Marqu ardt minimization procedure (iLevenbergI 
119441: lMarquardail968h is used extensively in many areas where 
the minimum of a nonlinear function needs to be located. It uses 
an adjustable mixture of a fully linearized solution and a steepest 
descent gradient search to locate the global minimum of a func- 
tion of an arbitrary number of variables. It's robust convergence 
properties make it the method of choice for many atmospheric 
inversion codes, and it is also at the core of the SPINOR code. 

The method is built on a linearization of the fit-quantity (the 
spectrum) around a "best guess" set of fit parameters (the atmo- 
spheric parameterization parameters). It requires computing the 
derivative of each synthesized data point yi with respect to each 
fit parameter of/, each of which composes one element of the so- 
called Jacobian matrix J. 

Although the calculation of J is generally possible 
only by using a finite-diff'erence approximation with re- 
spect to the fit quantity, the specific form of radiative 
transfer equation (RTF) allows us to compute them by 
moving the derivatives inside the formal soluti on inte- 
gral (iLandi Degl'Innocenti & Landi Degl'Innocentil 1 19771: 



iRees et al.l Il989l: iRuiz Cobo & del Toro Iniestal Il992h and 
evaluating the RTF using the derivatives of the opacity and 
emissivity with respect to the atmospheric quantities instead of 
the opacity and emissivity themselves. The resulting response 
functions (RFs) can be computed at a relatively modest addi- 
tional cost alongside the emergent spectrum and correspond 
to the partial derivatives of the spectrum with respect to the fit 
quantities, required by the Levenberg-Marquardt minimization 
procedure. 

Assuming the response functions are known, the solution 
a = must describe the observed data d 

Ja = d. 

For an approximate solution a' 

Ja' = d 9^ d. 

so that by using the linearity of J we may write 

J(a -a') = Jy = d-d'=p, (2) 

where v is the linear correction for the set of atmospheric param- 
eterization values a and JS is the diff'erence between the observed 
data and the spectrum synthesized using a. Since we generally 
have more data points then fit parameters, our desire to invert J 
directly to obtain y is frustrated by J not being square. The so- 
lution is to multiply both sides by the transpose matrix J^, thus 
creating the pseudo-inverse problem 

fjy = (3) 

where multiplying p by can be seen as converting the overde- 
termined but exact problem into a least square fitting problem. 
Although this linear system gives us the least-square solution 
to the linearized problem, unfortunately, the original inversion 
problem is usually far from linear, so that direct application of 
this solution may be anything but an improvement. An additional 
damping constant Am is therefore added to ([3]) to control the be- 
havior. The result is the familiar Levenberg-Marquardt form 

(J^J + ^Mdiag[J^J])y = S, (4) 

with S = J^J3 and Am the Marquardt damping constant. The ma- 
trix on the RHS can now be inverted to solve for the correction 
vector y. It is easy to see that for increasing values of the damp- 
ing constant Am, the fit matrix becomes more and more diago- 
nal, while the corrections become smaller and smaller, thus ap- 
proaching the steepest descent minimum search, which is known 
to be very slow and very vulnerable to finding a local minimum 
that is not the global minimum, but is guaranteed to lead to a 
decrease in the value of the merit function. 

The convergence strategy is now clearly to keep the damp- 
ing parameter as small as possible, in order to avoid the slow 
convergence properties of the steepest descent method, but large 
enough not to overstep the minimum so far as to increase the 
merit function. Ideally, the initial steps should only be concerned 
with the large-scale behavior of the minimization function; 
small-scale details should only come into focus once the domi- 
nant features have already been fairly well fitted. Unfortunately, 
as with all blind search algorithms, no guarantee exists that the 
minimum that is located is the lowest minimum. 
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2.3. Spatial coupling 

Now that we understand the way in which the Levenberg- 
Marquardt method optimizes the merit function, we are ready 
to consider the spatial coupHng between the pixels in an actual 
dataset. 

There are at least two important ways in which the spectra in 
adjacent pixels depend directly on each other: horizontal radia- 
tive transfer effects and instrumental scattering and diffraction 
effects. The first involves the considerable task of solving the 
RTE in non-LTE and in full 3-D and computing the response 
functions for all image elements and all fit parameters in these 
image elements. Although it is not clear a priori that this ef- 
fect does not have to be consider in general, we assume here 
that there are at least some photospheric lines for which the as- 
sumption of LTE is approximately valid so this effect can be 
neglected. The second way is easily seen to constitute a much 
simpler problem than the first, since the local derivative of the 
spectrum to the atmospheric fit parameters is simply spread to 
the neighboring pixels according to the PSF. 

Taking a closer look at the Jacobian matrix of a spectro- 
polarimetric image makes it clear that for spatially uncoupled 
inversions, without a PSF, the Jacobian is simply an assembly of 
uncoupled sub-Jacobians: 



r J 



1,1 



^ 



'1,2 



m-l 





so that we simply have 

f V 

^1,1 



(5) 











" n n 



(6) 



and the required fit matrix J^J clearly remains block-diagonal 
and can be solved efficiently block-by-block. 

For spatially coupled inversions with a uniform PSF (f(x, y) 
that only depends on the spatial coordinates, the Jacobian can be 
derived by making use of }T]) and can still be written in a blocked 
form: 



^o,oJi,i 

<^0,lJl,l 



^0-lJl,2 
<j^O,oJl,2 



<Pl-n,2-mJn,m-l <Pl-n,l-mJn,m 
^l-n,3-mJn,m-l ^l-n,2-mJn,m 



<^n-l,m-2Jl,l <^n-l,m-3Jl,2 
V^n-l,m-lJl,l ^n-l,m-2Jl,2 



^0,0 J n,m-\ 
^0,lJn,m-l 



^0,-1 J n,m 



but is clearly no longer block-diagonal. In practice, the extent 
over which (p does not vanish is much smaller than the com- 
putational domain, resulting in a relatively sparse system. The 
transpose is readily written down as 



<^o,oJ[,i 

<^0,-lJ 



1,2 



<^0,lJ[i 
<^0,oJi^2 



^n-\,m-3^\2 ^n-\,m-2^\2 



^l-n,2-mJfifn-l ^l-n,3-7«J„ ^_ 

JT tT 
. ^ ... n,m ^l-n,2-mJn,m 



^o,ojr,^_i 

^0,-1 Jj,m 



which, using (f (f = Y, allows us to write the desired matrix as 



5^0,0 Jfi J 1,1 

io,iJ^^2Ji'i 



^l-«,l-m J I i^i 



I i*3n,m 

T 



Yn-l,m- 
^n-l,m- 



Jl 

1,1 



'n^m-l^n,m 
^0,0 J« m^n,m 



^0,-1 



(7) 



which contains the matrix product of all jjj with all Jx,y, 
weighted with the autocorrelation function of the PSF centered 
on the location of the J;^,^;. 

If the PSF is not uniform, it is straightforward to replace the 
coefficients (pij with appropriate vectors containing the elements 
of the appropriate PSF. Although the computational overhead of 
this is probably not prohibitive, we leave this to be developed in 
a future paper. 

The Herculean effort of directly inverting this linear sys- 
tem, containing (Ua Xris XrixX Uyf elements, must not only be 
considered well out of reach of any computing resource likely 
to be available in the near future, if it were accomplished by 
some means, it is very likely that the large number of operations 
needed for computing the solution vector would lead to acumu- 
lation of unacceptable numerical errors, since the inverse matrix 
would likely not be sparse. 

However, since the forward problem is sparse, we are able to 
accurately evaluate 



(J^J - ^diag[J^J])x, = AXn = dn 



(8) 



and compute the defect d„ in the current estimate x„ of the cor- 
rection vector y 

d. =^-^„ (9) 

in a reasonable amount of time. We now approximate A with a 
sufficiently accurate approximation. A*, which is easy to invert 
but captures the essence of A. The solution of 

A*^n = dn (10) 

can now be used to improve the estimate 

repeatedly, until | d„ p< 6, with e a suitably small, positive num- 
ber. This process converges rapidly if an accurate approximation 
of A is used. 

2.4. Approximate operator 

In the process outlined above, an approximation to the full ma- 
trix A = J^J - /lMdiag[J^J] is required that captures the basic 
properties of the full operator but is both cheap to invert and has 
a sparse inverse. The block diagonal approximation is an obvious 
candidate, because it works well for a compact PSF, but requires 
a large number of iterations to converge when a more extended 
PSF is used. 

An extended-block variant of this method was also tried, 
where a spatially connected region was selected and inverted 
explicitly, as illustrated in Figure [T] This method generally im- 
proves the convergence characteristics considerably, but comes 
at the increased cost of inverting a larger part of A explicitly. 
Apart from the increased cost of the inversion, this approach suf- 
fers from the limitation that when the size of the blocks becomes 
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Fig. 1. Full matrix for a grid of 6x6 pixels (left), approximated 
by 9 blocks of 2 x 2 pixels (middle) and 4 blocks of 3 x 3 pixels 
(right). The cost increases rapidly with the number of pixels N, 
as the inversion of the resulting matrix scales oc N^.ln the upper 
right corner, the spatial distribution of the pixels is indicated by 
the color of the block they belong to. The full matrix has in ad- 
dition the cross-correlation function, Y, inserted in the lower left 
corner. 

too large, the numerical errors in the inverted block matrix can 
become large enough for the iterative approximation to fail to 
converge completely. It is therefore important to find the opti- 
mum balance between the cost of explicit inversion and the cost 
of applying a larger number of iterations, which generally de- 
pends on the problem. 

Other variations on the extended block method can be made, 
when more nonlocal behavior needs to be captured adequately. 
The staggered block method from Figure [2l can be useful when 



\\\ ■ 

















Fig. 2. Full matrix for a grid of 8x8 pixels (left), approximated 
using a staggered configuration of 2 x 2 pixel blocks (right). The 
approximated matrix has more than 9 blocks, due to boundary 
eff'ects, compared with only 9 above, but many of the blocks are 
smaller. The insets are as in Figure [T] 



interaction over long distances is important, such as for sig- 
nificantly oversampled data. Intermittent application of several 
methods is possible, but the need to recompute the inverse of 
the approximate matrix generally outweighs the possible gain in 
convergence. 

Computational optimization using explicit storage of all 
nonzero elements of (J^J-Am diag[J^ J]) and the approximation 
of the inverse, in combination with multithreading of the itera- 
tive solver, reduces the additional time overhead of this method 
to only ^^30% of the cost of synthesizing the Stokes profiles and 
response functions with the SPINOR code, for critically sam- 
pled data and with observational data containing -100 spectral 
points. This overhead is expected to decrease further when the 
amount of spectral points increases. 

Depending on the approximate operator that is used, the con- 
vergence of the inversion of A typically requires on the order of 
100 iterations to reach an error in the defect of 10%, as shown 



in Figure [3l This takes significantly less time than the computa- 
tion of A and A* and so does not significantly contribute to the 
overall execution time, as long as the inverted submatrices can 
be stored in RAM. 

Unfortunately, the need to update the correction vector and 
recompute the defect, requires knowledge of the current estimate 
across the extent of the Y functions, resulting in the need for a 
significant amount of communication when a distributed paral- 
lelization strategy is desired. For this reason, at present, only an 
efficient multithreaded algorithm exists, requiring a shared mem- 
ory environment. 

In some cases, additional shortcuts can be taken to reduce the 
work requirements. Such shortcuts, however, do not always re- 
sult in unconditional convergence, as is the case for the complete 
treatment as outlined above. One particularly useful approxima- 
tion in the case of a PSF with very extended wings is to calculate 
the defect J3 explicitly using the complete PSF, while calculating 
the Levenberg-Marquardt correction y using only the central part 
of the PSF. This results in good convergence when the value of 
the PSF integrated over the extended part of the PSF does not 
exceed the value of the PSF integrated across the central part. If 
this condition is not met, an exponentially growing correction is 
obtained and the inversion process will diverge. 

2.5. Inversion strategy 

Although we can evaluate the global problem, it is not self- 
evident that this is the optimal way to proceed. If we consider 
the pixel-by-pixel approach as completely uncoupled and a very 
extended PSF as a global inversion problem, it is clear that in 
general all levels of coupling will be encountered. For the calcu- 
lations presented in this paper, only the global problem has been 
considered; however, it should be kept in mind that the global 
approach may not be the optimal choice in all cases. 

The value of the merit function, is plotted in Figure [3] for 
the inversion of an MHD simulation as a function of the iteration 




20 40 60 80 100 120 140 10° 10^ 10^ 

Iteration Iteration 

Fig. 3. Convergence properties of the iterative step determination 
and the global inversion left: error in the defect as a function of 
iteration number, right: average (black) and maximum (red) 
as a function of iteration number. 

number. The convergence rate can be strongly influenced by em- 
ploying additional smoothing at regular intervals, which can ef- 
fectively push the solution out of local minima by damping out 
large local differences that typically result from the cancellation 
of errors by other errors with the opposite sign. 

A returning issue when converging oversampled data is the 
development of a large number of local minima, frequently in- 
volving spatial frequencies in the solution exceeding the diffrac- 
tion limit. One of the more effective ways to reduce the proba- 
bility of getting stuck in one of these minima is to initially limit 
the corrections applied to the solution to the lower frequencies 
and to relax this limitation as the inversion progresses. 
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Another effective way to avoid a slow and cumbersome de- 
scent to the global minimum is to introduce a significant amount 
of noise in the data in the first stages of the inversion. The noise 
eff'ectively hides the less significant spectral features from the 
view of the inversion algorithm, leading to a rapid first descent 
phase. After the minimum has been approached with sufficient 
accuracy, the noise is gradually reduced, allowing for the details 
to enter the solution. This method is fairly effective in avoiding 
getting captured by outlying local minima, although it is not very 
effective in avoiding false minima near the global minimum. 

3. Simulations 

To evaluate and optimize the performance of the spatially cou- 
pled inversion problem, it is necessary to know what the value 
of a particular atmospheric parameter really is, something that 
can only be done in the context of simulated atmospheres. The 
atmosphere chosen for the tests below is a 288x288x100 point 
MURaM MHD simulation dVogler et al . 2005) of the quiet Sun 
with a mean field of around lOOG and a spatial resolution of 
20.8km (0.029"). The grid size allows for convenient binning in 
a number of coarser grid sizes, which we make use of in the tests 
carried out below. 




X [Mm] 



Fig. 5. The continuum intensity of the employed MURaM quiet 
Sun simulation. 

The verification procedure necessarily remains strictly lim- 
ited to validation and verification of the method under inherently 
ideal circumstances, since the forward and the backward prob- 
lems use the same physics, a situation that is not likely to oc- 
cur in the case of real observations. The effects of errors in the 
physics will be considered inherent to all inversions and are left 
for a later discussion. 

The starting point of our calculations is an MHD cube, pic- 
tured in Figure \5\ that provides the relevant atmospheric quan- 
tities (temperature, pressure, velocity and magnetic field) as a 
function of three spatial coordinates. An unambiguous compar- 
ison of such a fully stratified atmosphere with an inverted sim- 
plified atmosphere is a complicated task by itself and one that is 
not unique to the coupled inversion problem. 

To steer clear of a complicated discussion on this topic, the 
atmosphere was first used to synthesize the emergent spectra, 
which were then inverted using simple three node atmospheres, 
with the nodes at the same heights as they will be for the in- 
versions. The resulting atmospheres are only approximations to 
the original ones, but they have the advantage that they can be 
represented exactly by the inversion code. 

The fitted atmospheres are then used to generate a second 
data cube of emergent spectra, with the notable difference that 



these spectra can be reproduced exactly by the simplified three- 
node atmospheres used by the inversion code, so that it should be 
possible to invert the spectra perfectly to the input atmospheres 
used to generate them, provided the mapping of the atmosphere 
to the emergent spectra is unique. 

3.1. Ideal atmospheres 

As motivated and explained above, we attempt to validate the 
method under the simplest possible circumstances, where the 
profiles are generated from a known solution that can be repro- 
duced exactly by the inversion code. To achieve this, the three- 
node atmospheres, fitted to the synthesized spectra as described 
above, were first binned to the resolution of the artificial dataset 
we want to generate, before the spectra to be degraded and in- 
verted are synthesized from them. 

If the inversion problem is well-posed, this simplification 
should allow the inversion code to recover the binned atmo- 
sphere with arbitrary precision in the absence of spectral and 
spatial degradation. Since this is generally not the case for all 
pixels in the field of view (a few spectra will always be diffi- 
cult to invert), to evaluate the properties of the coupled inversion 
as independently of the properties of the inversion procedure as 
possible, we focus on the differences between the pixel-by-pixel 
inversion of the undegraded data and the coupled inversion of 
degraded data. 

3.1.1. Pixel-by-pixel inversions 

After a first verification that the atmosphere used to generate the 
spectra gives a perfect match to the artificial data, a pixel-by- 
pixel inversion was carried out of the whole FOV, to test the in- 
herent invertibility of the atmospheres using the algorithm used 
by the inversion code. Interestingly, the results indicate that, al- 
though the code should be able to perfectly reproduce the spec- 
tra, it does not succeed in doing so for approximately 1 % of the 
pixels in the FOV. 

Repeated "nudges" of the atmospheric parameters in random 
directions in parameter space improve the overall success rate, 
but in a number of pixels, the solution is evidently difficult to 
locate using the minimization procedure employed by the code, 
so that significant differences between the artificial observation 
and the inverted spectrum remain. However, where the inversion 
code indicates a good fit, the atmospheric parameters are accu- 
rately retrieved, suggesting that the inversion problem is not un- 
derdetermined in this simplified situation. 

The fraction of the pixels that does not produce an accurate 
inverse is approximately independent of the resolution, which is 
somewhat surprising, since although the atmospheres are binned 
to the resolution of the data, so that no unresolved structure can 
be present, the profiles tend to exhibit more extreme behavior at 
higher resolution. Figure |4] shows the retrieved inverted values 
as a function of the input values. The temperature and Doppler 
velocity appear to be very well-determined in all layers of the 
atmosphere, but the magnetic field strength and orientation show 
minor scatter around the diagonal, indicating that there are limits 
to the invertibility of the problem. 

Figure [6] shows the recovered Stokes parameters as a func- 
tion of the input Stokes parameter. The agreement is very good in 
all Stokes parameters, indicating that the minor scatter observed 
in the recovered magnetic field parameters is due to [near] de- 
generate behavior in these parameters, possibly owing to an in- 
herently low response. 
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Fig. 4. 2-D histograms of inverted atmospheric parameters (vertical) as a function of input atmospheric parameters (horizontal), for 
a coupled inversion of simulated degraded 50cm telescope data (top row) and an uncoupled inversion of identical but undegraded 
simulated data (bottom row), for 5 atmospheric parameters across all heights in the atmosphere. 
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Fig. 6. 2-D histograms of inverted Stokes parameters (vertical) as a function of input Stokes parameters (horizontal), for a 2-D 
inversion of degraded 50cm telescope data (top row) and a 1-D inversion of identical but undegraded data (bottom row). 



3.1.2. Spatially coupled inversions 

In the spatially coupled case, we must choose a PSF to degrade 
the spectral data with. With a large volume of available po- 
tential da ta in mind, we choo se the central part of the PSF of 
the SOT (Tsiineta et al.ll2008l) on board the Hinode spacecraft 
(iKosugi et al.ll2007h . as indicated in Figure [TOl so we get a direct 
indication of the performance we can expect for such data. 

The original resolution of the MHD cube, 0.029"/pixel, is ap- 
proximately a factor 9 higher than the diffraction limit of 0.26" 
of the Hinode SOT at 6300A, so that a binning factor of 5 would 
be the closest approximation of critical sampling. Since this is 



not an integer factor of 288, the cube was binned by a fac- 
tor 6 to 48 X 48 pixels, corresponding to 0.174"/pixel, which 
is slightly undersampled and in r easonable agreem ent with the 
Hinode spectro-polarimeter (SP) (iLites et al .1120011) pixel size of 
0.16". In addition, a cube binned by a factor 3 to 96 x 96 pix- 
els was produced to investigate the response of the algorithm to 
oversampled data with unresolved fine structure. 

Compared to a pixel-by-pixel inversion, the spatially cou- 
pled inversions clearly require an increased amount of human 
intervention, since there are so many coupled fit parameters that 
there is a significantly increased probability of finding a local 
minimum that is not the global minimum. Identifying problem- 
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Fig. 7. Visual comparison of a selection of inverted parameters from spatially degraded spectra with the original data (top) for a 
spatially coupled inversion (middle) and a spatially uncoupled inversion (bottom). Shown are the parameter maps at the optical 
depths where they are most accurately determined. From left to right: Temperature at r = 1 , magnetic field strength, inclination 
angle and azimuth angle at r = 0.16, and the line-of-sight velocity at r = 1. The color schemes are optimized for each quantity 
and have an identical but otherwise arbitrary scale for the original and inverted parameters. Tick marks indicate the spatial scale in 
arcseconds. Clearly visible is the diff'erence in the accuracy of the retrieved parameters when the spatial smearing by the telescope 
is not taken into account. 



atic pixels is difl&cult, because the only measure of the quality 
of the fit, the diff'erence between the convolved profiles and the 
data, is the result of an ensemble of profiles, so that it is not pos- 
sible to rate the contribution to the fit quality of one particular 
profile. A "nudge" of all fit parameters in a random direction is 
not necessarily helpful, since many pixels may have already con- 
verged and a random step in that solution may lead to improved 
convergence in one place, but not in another. 

For slightly undersampled data, the method does not seem 
to have much more trouble recovering the solution than in the 
case of a 1-D inversion in the absence of any spatial degradation 
for the majority of the atmospheric parameters. Figure |4] shows 
the retrieved atmospheric quantities as a function of the input 
atmospheric quantities used to generate the undegraded spectra. 
Although the scatter around the diagonal is somewhat larger than 
for the pixel-by-pixel inversions of spatially undegraded data, 
when we keep in mind that the histograms in Figure |4] have a 
logarithmic gray scale, we may conclude that on a linear scale, 
the only parameters with significantly increased scatter are the 
horizontal components of the magnetic field. 

This is also clearly visible in the inverted Stokes spectra 
shown in Figure [6j where the recovered Stokes Q and U val- 
ues show the greatest deviation from their real value. Since the 
value of the merit function does not vanish for this solution, we 
can state with confidence that despite the human intervention. 



the inversion code found a local minimum that is not the global 
minimum, which is the most likely cause for the observed errors. 

As a final illustration, Figure|7] shows the maps of a selection 
of the original parameters used to calculate the Stokes profiles, 
together with the inverted result. To emphasize the importance 
of correctly treating the spatial degradation, the uncoupled (ID) 
inversion result of the spatially degraded data is shown. Even 
though the parameters are shown at the optical depth where they 
are most accurately constrained by the data, the loss of accuracy 
caused by neglecting the spatial smearing of the data is clearly 
visible in the inverted results. We will explore the diff'erences 
between the various inversion strategies in more detail in a future 
paper specifically on this topic. 

3.2. The effect of noise 

The eff'ect of noise is one of the main problems to deal with when 
one is interested in the interpretation of real data. To investigate 
the eff'ect of additive Gaussian noise on the spatially coupled in- 
versions, a number of synthetic observations were produced with 
noise levels that are typical for data obtained with an efficient in- 
strument and an integration time of several seconds. 

The results of a repeat of the test from section 13.1.21 using 
data with varying levels of noise, is shown in Figure[8j The criti- 
cal level for reliable recovery of the atmospheric parameters ap- 
pears to be 10"^ Ic, beyond which first the transversal (T) mag- 
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Fig. 8. Two-dimensional histograms of inverted atmospheric parameters (vertical) as a function of input atmospheric parameters 
(horizontal), for a 2-D inversion of degraded 50cm telescope data for noise levels of (left to right) 0, 10"^, 10"^, 3 x 10"^, and 10"^ 
/c, for 5 atmospheric parameters across all heights in the atmosphere. 



netic field properties rapidly degrade, followed by the longitudi- 
nal magnetic field and finally the temperature and line-of- sight 
(LOS) velocity. At a noise level of 10"^ /c, practically all at- 
mospheric parameters show large errors, suggesting that a level 
well below that should always be achieved for any atmospheric 
property to be reliably recovered. 

In particular, the recovered temperature and line-of- sight ve- 
locity appear to only be weakly aff'ected by the noise. This is 
easily understood as these quantities may be inferred directly 
from the intensity, which is clearly the Stokes parameter least 
aff'ected by the noise. In addition, the intensity in the continuum 



has a strong response to the temperature in the deepest layers, as 
well as a large number of spectral data points constraining that 
part of the profile. 

The atmospheric quantities that only depend on the Stokes Q, 
U, and V parameters are significantly more sensitive to the noise, 
which is not only much more severe owing to the relatively low 
signal level of these data, but also to the negligible response of 
the majority of the (continuum) data points to these quantities. 
The final consequence of this is that for reliable recovery of the 
longitudinal component of the magnetic field in the quiet sun, a 
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signal-to-noise of ~ 10 ^ Ic is required, for the orientation and 
field strength it must be significantly better. 



3.3. Unresolved structure 

One of the main concerns in any inversion is the uniqueness 
of the solution, in other words, is the problem invertible at all? 
Naively one might expect the spatially coupled fit to have prob- 
lems here, since only the sum of a number of spectra is con- 
strained, not each individual spectrum. The difficulty the code 
has in locating the global minimum solution suggests, however, 
that this is actually not the case, so that if there are multiple op- 
timum solutions, clearly none of them is found. 

In an attempt to test what happens to the solution in cases 
where this uniqueness should break down, an inversion was done 
on oversampled simulated data, a case where conventional con- 
volution is known to no longer provide one-to-one mapping. To 
this end, a data cube was simulated as before, but this time at a 
resolution of 0.085"/pixel, and convolved with the Hinode PSF 
sampled at 0.085"/pixel, so that it is possible for the solution to 
contain information that is not constrained since it is no longer 
present in the inverted observation after convolution with the 
PSF. 

Not unexpectedly, the inversion does not find as good a so- 
lution as in the case of critically sampled data, even if just the 
value of the merit function is considered, but it shows no sign of 
instability or amplification of errors either. This suggests that the 
constraints provided by imposing a particular atmosphere and 
particular physics on the solution are able to stabilize the inver- 
sion sufficiently, even in oversampled cases. 

Figure [9] shows the scatter of the inversion result as a func- 
tion of the input atmosphere. Although the scatter is significantly 
increased over the critically sampled case, it is similar to the 
scatter of an appropriately spatially filtered atmosphere, shown 
in the bottom half of the same figure. The algorithm apparently 
does not recover much information beyond the diffraction limit, 
but instead appears to recover an atmosphere that has all fre- 
quencies beyond the diffraction limit removed from it. 



3.4. Discussion 

The simulation results show that the evaluation of the merit func- 
tion in convolved data space allows for accurate discrimination 
between profiles in the undegraded image to spatial scales up to 
the diffraction limit of the PSF, if the instrumental properties are 
known. The undegraded Stokes profiles are recovered with a mi- 
nor increase in the absolute error, a behavior closely followed by 
the inverted atmospheric parameters. Clearly, the quantities that 
depend primarily on the Stokes I profiles are recovered more ac- 
curately than those depending on Stokes Q, U, and V, even if 
only numerical noise is present. 

Interestingly, the inversion algorithm does not converge to 
the absolute minimum, where the merit function was verified to 
vanish. The difficulty is that whereas a pixel-by-pixel inversion 
allows the identification of poorly converged pixels, this is not 
possible in the spatially coupled case. 

The failure to find the absolute minimum within the space 
accessible by the algorithm suggests that if the solution is de- 
generate, the alternative degenerate solutions are as difficult to 
locate as the true solution, even in the case of oversampled data. 



4. Application to Hinode observations 

Since the spectral dimension is crucial for this method to work 
and the simplification introduced by the assumption of a con- 
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Fig. 10. left: The Hinode PSF at 0.16" per pixel resolution on a 
logarithmic gray scale from 10"^ (black) to 1 (white). The black 
box is the explicit part used in the spatially coupled inversions, 
right: Hinode SOT pupil function. 

slant PSF is advantageous in keeping the computations tractable, 
the spectro-polarimetric scans of the Fe I lines at 6301. 5A and 
6302. 5 A provided by the Hinode spectro-polarimeter were used 
to test the performance of the coupled inversion method on real 
solar data. 

The significant scanning time of this instrument compared 
to the evolution time scale of some solar structures gives reason 
to be critical of the result, however, the slight under- sampling of 
the data keeps the evolution effects of the solar image relatively 
small on the critical spatial scale of the PSF. Nonetheless, it is 
important to remember that accurate reproduction of the spectra 
may not be achieved in some cases. 

In addition to this obvious shortcoming, we are dealing with 
a real solar atmospheric structure, with, more likely than not, 
a complicated depth stratification, which will not be accurately 
represented by the simple atmosphere imposed by the inversion 
code. Furthermore, although presumably very accurate, the as- 
sumption of LTE probably does not hold exactly and the atomic 
and molecular line data are not known perfectly. 

We nonetheless proceed by computing the PSF of the tele- 
scope. The PSF was computed from the Hinode pupil as s peci- 
fied i n detail in the technical reference manual (Suematsu et al.l 
20081). As the f ocus p osition was found to differ from by 
Danilovic et al. (l2008h . a 0.1 wave defocus was added, even 
though the focus position of this dataset was not accurately 
known. The PSF up to a radius of 4.8" is shown in FigurefTOlon a 
logarithmic grayscale from 10"^ to 1, showing a complicated and 
extended interference pattern with significantly enhanced radial 
structures produced by the triangular spider. 

To calculate a 2D inversion using the whole PSF, an unrea- 
sonably large amount of memory and computing time would 
be needed. Therefore, the shortcut described in section 12.41 was 
used, where the central part of the PSF explicitly used in the 
inversions is indicated by a black square in Figure [TOl Since the 
wings of the PSF are taken into account explicitly when calculat- 
ing the defect and the merit function, they are properly included 
in the converged result. 

Since the results did not really require it, no straylight con- 
tamination was included in the inversion, which is not to say that 
there is none in the data. However, since the fit to the Stokes pro- 
file with the lowest intensity in the inverted data cube, which was 
found to be around 0.12Ic,hsra in the darkest part of the umbra. 
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Fig. 9. Two-dimensional histograms of inverted atmospheric parameters (vertical) as a function of input atmospheric parameters 
(horizontal), for a 2-D inversion of degraded 50cm telescope data (top row) and a spatially filtered but otherwise undegraded 
atmosphere (bottom row), for 5 atmospheric parameters (horizontal) across all heights in the atmosphere. 



was still accurate to within ~ OAIcjocah the amount of stray- 
light present that cannot be explained with the telescope PSF is 
unlikely to be more than around 0. 01 1 c,hsra, although no eff'ort 
was made to more precisely determine this number. 

A scan of AR10933, recorded on 5 January 2007, 11:54- 
12:27 UT was selected, showing a sunspot umbra, penumbral 
filaments, and some solar granulation, very close to the cen- 
ter of the solar disk. The FOV was selected since it contains a 
wide variety of atmospheric conditions and fine structure that 
are likely to benefit from a spatially coupled inversion. As in the 
simulations from section[3l three nodes were placed at log r = 
-2.5, -0.9, 0.0 to represent the solar atmosphere. 

4.1. Inverted observations 

Before the inversion results are discussed, it should be stated that 
the focus of this paper is a discussion of the coupled inversion 
method and not a discussion of the atmospheric structure that 
is the actual result of the inversion. Since discussing the latter 
would draw our attention away from the method, we leave a de- 
tailed analysis of these results to a separate paper and instead 
limit ourselves to a brief, qualitative assessment of the results 
combined with a comparison of the fitted undegraded profiles to 
the observed (degraded) data. More specifically, we produce im- 
ages generated from the inverted profiles (inverted images) and 
compare them to images made from the observed data. 

The inverted images at 6302. 55A of the selected active re- 
gion are shown in Figure [TT] The synthetic images generated 
from the inverted atmospheric structure (right side) compare 
very favorably with the original data (left side), in that they show 
much more contrast and a significant increase in the level of de- 
tail. This is even clearer in Figure [121 where the upper righthand 
corner of the umbra is shown enlarged for both the observed and 
the inverted data. 

A selection of atmospheric parameters is shown in Figure[T3l 
where the level of detail suffices to see the temperature drop in 
the penumbral dark cores at the highest node, the clear drop in 
magnetic field strength in the penumbral filaments, the large in- 
crease in inclination in the penumbral filaments, and clear signs 



3 
2 
1 

CO 

1 ° 

> -1 

-2 



\ 




1 1 ^1 1 1 1 1 

1 1 1^1 1 1 1 


1 1 ' 1 1 1 1 1 
1 1 1^ r 1 1 1 _ 



-3-2-10 1 2 3 

X [arcsec] 



-3-2-10 1 2 3 

X [arcsec] 



Fig. 12. Observed (left) and inverted (right) Stokes images in 
Stokes I (top) at and Stokes V (bottom) at 6302.55A. 



of downflow on both sides of the filament in the deepest layers, 
only seen before in datasets taken with larger solar telescopes 
(iJoshi et al.ll20TTl: IScharmer et aDl201lh . 

Despite this obvious improvement in the inverted parame- 
ters, the results also show several regions with obvious signs of 
oscillations on the scale of a single pixel. The oscillations show 
up in specific regions and cannot be removed by offering the 
code a smooth start condition. The pattern itself is confined to a 
certain region and does not spread beyond that, but its phase is 
not fixed. Several inversions starting from slightly different start 
conditions all produce a similar oscillatory pattern in the same 
region, but none of them is identical to any of the others. An av- 
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Fig. 11. Stokes images in the Fe I 6302. 5A line+ 57 mA, before (left half) and after (right half) inversion for the Stokes parameters 
I (top left), Q(top right), U(bottom left) and V(bottom right). 



erage of several such solutions does not show any oscillations, 
but if an inversion is started from this averaged solution, it will 
develop new oscillations. 

The regions in which these oscillations are found are typi- 
cally fine structured and fast evolving, such as the outer penum- 
bra. One example is shown in Figure [141 where the oscillations 
are clearly visible in the undegraded Stokes V and the LOS ve- 
locity in the penumbra, but are perfectly smooth in the granu- 
lation right next to it. The degraded Stokes V image shows that 
these oscillations average out to form a smooth image with only 
very weak oscillations. 

The persistence of the oscillations in combination with the 
relative insensitivity of the inversion to their phase, but not their 
amplitude, suggest that we may not be dealing with just any de- 
generate solution here, but that the oscillations may actually be 



required to produce a spectral feature that cannot be produced 
by the atmospheric model. In particular, since the occurrence ap- 
pears to be concentrated in the deep layers, the code must have 
difficulty accurately reproducing the wings of the line profiles. 
This may be due to calibration errors, but it is more likely caused 
by fine-scale structures that are not resolved by the telescope, 
possibly in combination with significant time evolution effects 
occurring during the spatial scan. 

Apparently the inversion algorithm prefers to compensate for 
these errors by introducing what is in effect a microturbulent ve- 
locity structure, that is, by alternating the LOS velocity on the 
smallest possible scale. The averaging occurring in the convolu- 
tion of the synthesized profiles then produces the broadened pro- 
files that are observed. Interestingly, the removal of the micro- 
turbulent velocity as a fit parameter significantly increases the 
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Fig. 13. A selection of inverted atmospheric parameters. From 
left to right, top to bottom: T-2.5, ^-0.9, T-0.9, and Vd,o, with 
the subscript denoting the value of ^^logr. The color scale used 
for r_2.5, ^-0.9, and 7.0.9 is indicated on the top right, the one 
for Vd on the bottom right. The ranges of the color scale are 
[3300-5200] K, [500-2500] G, [0-180]°, and [-3,3] km s-\ re- 
spectively. 



problem, suggesting that this parameter is used to fit some of the 
broadening caused by unresolved structure, but that significant 
diff'erences remain that cannot be reproduced in this way. 
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Fig. 14. Inverted Stokes images in Stokes I and V (left), along- 
side with the inverted LOS velocity at r = 1 and the degraded 
Stokes V. Clearly visible are the oscillations in the penumbra, 
which are invisible after degrading the image. 



The temperature map in Figure [13] shows some horizontal 
stripes that can clearly not be associated with real solar features. 
They are parallel to the spectrograph scanning direction and can 
only be attributed to imperfections in the flatfields, made visible 
here by the level of detail to which the observed data are fitted. 
These artifacts are concentrated in the top layer of the inverted 
atmosphere, most likely owing to the small number of spectral 
points that are available to constrain this part of the atmosphere 
and the intrinsically lower intensity that these points have com- 
pared to the other parts of the spectrum. 

Even a small calibration error in the width of the instrumen- 
tal transmission profile will cause compensation eff'ects in the 
solution, predominantly in the top layer. Indeed, a deliberate un- 
derestimate of the spectral transmission profile of the instrument 
introduces significant oscillations in the top layer of the atmo- 
sphere, as the algorithm attempts to broaden the line by alternat- 
ing an atmospheric quantity. 

It is clear that the algorithm is able to successfully fit the ob- 
served spectra, but it is able to do so with such high accuracy 
that it is essential to calibrate the data to a precision level that 
matches or exceeds that of the inversion result, which is signifi- 
cantly higher than the noise level of the data. 

4.2. Deconvolved observations 
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Fig. 15. Deconvolved (left) and inverted (right) Stokes images in 
Stokes I (top) at 6302. 55A (line wing) and Stokes U (bottom) at 
6302. 65A (near continuum). Note the enhanced contrast in the 
high signal case (top) and the improved noise characteristics in 
the low signal case (bottom) of the inverted result as compared 
to the deconvolved result. 

To evaluate the results obtained with the coupled inversions, 
the inverted undegraded profiles were compared to the data de- 
convo lved u sing Richardson-Lucy deconvolution (iRichardsonI 
19721: lLucvl lT974). Figure [15] shows two images, one decon- 
volved and the other one the result of the coupled inversion. 
The Stokes I images, at 6302. 55A, show a significant diff'erence 
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Fig. 16. Angular averaged power spectrum of the deconvolved 
(black) and inverted (red) solutions of Stokes I at 6302. 55A, 
approximately half way into the Fe I line at 6302. 5A. Clearly 
visible is the increased power of the inverted result in the high- 
frequency range of the power spectrum. 



in contrast, as is confirmed by the azimuthally averaged power 
spectrum shown in Figure O indicating that the inverted im- 
age contains much more power at the high frequencies than the 
deconvolved result. 

A low-signal image, on the other hand, such as the Stokes U 
image shown at the bottom of Figure \T5\ shows a power spec- 
trum that has much more power than the inverted result, as can 
be seen in the righthand panel of Figure [161 Close inspection of 
the deconvolved image confirms that this can be attributed to a 
high level of noise in the image, not to any actual structure. The 
inverted result, on the other hand, still shows credible, coher- 
ent structures with nine orders of magnitude less power than the 
maximum power in the Stokes I image and does not contain any 
noise. 

The inversion is able to produce high-frequency informa- 
tion when it is supported by the data, but suppresses it very 
eff'ectively elsewhere. The likely reason for this lies in the fact 
that the inversion is based on the complete data cube, whereas 
the deconvolved result only contains information from a single 
wavelength. The eff'ective noise level of the deconvolved result 
is therefore much higher than for the inverted result, specifically 
at very low signal levels, where the inverted result is constrained 
by the high- signal regions of the spectra, where noise is much 
less important. 

5. Conclusions 

We have developed a new method for inverting 2-D maps of 
spectro-polarimetric data that have been degraded spatially in 
a known way. The method uses the information contained in the 
spectral dimension and the known spatial degradation properties 
to constrain a parameterization of the atmosphere over the whole 
FOV simultaneously. 

The method was formulated in a manner specific to an inver- 
sion code using a "greedy" optimization method, and although it 
is implemented for the SPINOR inversion code, it should not be 
considered as specific to this code. The method only changes the 
way in which the merit function is evaluated and modifies the 
response functions accordingly, implying that it should be easy 
to adapt to other inversion codes. The suitability of the specific 
implementation described in this paper is limited to the specific 
category of minimum search algorithms and cannot be used di- 
rectly in Monte-Carlo or genetic codes without adaptations. 

Tests carried out on simulated MHD cubes show that on a 
critically sampled grid, the code is able to retrieve atmospheric 
parameters by inverting profiles synthesized from them and de- 
graded using a realistic PSF, with only a small increase in the 
error of the retrieved parameters as compared to a pixel-by-pixel 



inversion of the undegraded spectra. For oversampled data, the 
code is able to successfully recover the low-frequency compo- 
nents of the solution, but appears to have difficulty recovering 
information beyond the diffraction limit. At no time was unsta- 
ble behavior resulting from oversampling the solution observed. 

The method developed here is a specific case of a more gen- 
eralized approach to data reduction, involving the application of 
calibration information to synthetic data (e.g. in the forward di- 
rection), instead of "correcting" the real data with the inverse 
calibration. Thanks to this forward application of the calibration 
information, the noise amplification typically associated with the 
correction of degraded data is strongly reduced, since the in- 
verted calibration data are never applied to it. In addition, by 
combining all the available information to constrain a simplified 
parameterized atmosphere, all errors and limitations are accu- 
mulated in one place, resulting in a robust solution with consid- 
erably improved specifications over the input data. 

Based on the results of section |4l we conclude that this 
method is suitable for application to spectral imaging data, 
where the assumption of a constant PSF can be made. There 
is no fundamental problem extending this method to the treat- 
ment of data where this assumption is not valid, which will be 
explored in a future paper. 
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